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ABSTRACT. 

Numerical  simulations  of  the  weak  beam-plasma  instability  were  done  in 
the  turbulent  regime  where  small-scale  trapping  is  a  dominant  feature  of  the 
instability,  a  regime  with  behavior  not  predicted  by  quasilinear  theory.  It  occurs 
when  the  trapping  frequency  u  =  ( k2D)^ is  larger  than  the  growth  rate  7* 
of  the  instability.  The  results  of  the  simulations  were  compared  with  those  of 
a  specific  model  of  the  turbulence,  the  so-called  "turbulent  trapping”  model, 
which  gives  precise  formulas  for  the  particle  correlation  functions,  and  predicts 
a  growth  rate  well  enhanced  over  the  quasilinear  value.  It  was  found  that 
the  model  gives  accurate  predictions  for  the  correlation  functions,  and  thus 
provides  a  good  description  of  the  turbulent  structure  of  phase  space.  On 
the  other  hand,  while  growth  rates  were  enhanced  over  the  quasilinear  values, 
the  enhancements  observed  are  smaller  than  expected  from  the  quantitative 
predictions  of  the  model.  Further  work  is  necessary  to  determine  whether  this 
discrepancy  is  a  failing  of  the  turbulent  trapping  model,  or  the  result  of  the 
numerical  limitations  of  our  computational  scheme. 
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I.  INTRODUCTION, 

Since  its  creation  over  20  years  ago  (*),  quasilinear  theory  has  been  ac¬ 
cepted  as  the  valid  description  of  weakly  turbulent  plasmas  where  wave-particle 
interactions  dominate,  with  as  central  paradigm  the  weak,  warm-beam-plasma 
instability.  It  is  only  relatively  recently  that  the  rigor  of  its  derivation  and 
its  range  of  validity  in  the  one-dimensional  case  have  been  questioned  by  J.C. 
Adam,  G.  Laval  and  D  Pesme  (ALP(,2,3,4,6)).  Their  work  demonstrates  the 
existence  of  a  new  regime  existing  within  the  accepted  limits  of  quasilinear  the¬ 
ory,  in  which  both  the  growth-rate  and  the  velocity-space  diffusion  should  be 
enhanced  over  their  quasilinear  values. 

This  conclusion  is  based  on  the  existence  of  strong  mode-coupling  effects 
between  resonant  waves,  due  to  the  self-consistency  of  the  electric  field,  ef¬ 
fects  which  occur  when  the  resonance  broadening  frequency,  ?/*,  s  (fc2!}9*)1/3, 
is  larger  than  the  quasilinear  growth-rate  7^.  In  this  regime,  a  set  of  mode 
coupling-terms,  arranged  in  an  infinite  series,  becomes  of  the  same  order  as 
the  quasilinear  term.  This  arises  because  of  the  partial  trapping  of  the  beam 
particles  by  the  waves  on  a  time  scale  (7fc  a  process  which  gen¬ 

erates  a  harmonic  series  of  sidebands  or  19 quasimodes”,  that  is  non-resonant 
perturbations  of  the  beam  distribution  function,  at  the  harmonic  frequencies 
uj  —  nwk  and  at  the  wave  numbers  nk.  These  quasimodes  do  not  satisfy  Pois¬ 
son’s  equation,  but  perturb  particles  for  which  the  resonance  condition  of  the 
fundamental,  u ;*  =  Art;,  continues  to  be  satisfied.  Because  of  this  invariance  of 
the  resonance  condition,  the  mode-coupling  coefficients  are  large,  even  though 
the  coupling  takes  place  via  small,  non-resonant  sidebands.  Thus,  the  cascade  of 
the  quasimodes  beating  back  into  the  fundamental  leads  over  a  few  e- foldings, 
to  effects  of  the  same  order  as  the  usual  wave-particle  interaction.  In  particular, 
under  this  regime  the  statistical  properties  of  the  electric  field  strongly  deviate 
from  Gaussian  statistics  and  this  results  in  turn  in  an  enhancement  of  the  dif¬ 
fusion  coefficient  over  its  quasilinear  prediction.  Through  energy  conservation, 
a  concommitant  increase  in  the  growth  rate  must  result  as  well. 

In  a  preceding  paper (6),  two  of  the  authors  presented  model  equations  that 
were  proposed  in  an  attempt  to  take  into  account  the  mode-coupling  terms 
which  originate  in  the  self-consistency  of  the  electric  field,  and  which  invalidate 
the  quasilinear  theory  in  one  dimension.  This  model  has  received  the  name  of 
the  "turbulent  trapping  model",  since  it  is  devoted  to  describing  those  physical 
effects  associated  with  the  partial  trapping  of  the  particles  in  the  wave  packets. 

In  the  present  paper,  we  numerically  study  the  interaction  of  a  very  weak 
warm  beam  with  a  cold  and  massive  bulk  plasma,  with  the  aim  of  detecting 
the  new  effects  predicted  by  the  turbulent  trapping  model.  Our  numerical 
experiments  are  simplified  as  much  as  possible,  but  hopefully  retain  all  the  basic 
physics  of  the  process.  The  most  fundamental  simplification  lies  in  modelling 
the  bulk  plasma  as  a  cold,  linear  fluid.  In  effect,  it  provides  nothing  more 
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than  the  linear  dielectric  which  sustains  the  electron  plasma  waves.  Apart 
from  a  small  numerical  effect,  these  are  waves  with  zero  group  velocity,  that  is 
oscillations  with  stationary  envelopes,  growing  on  the  energy  of  the  circulating 
beam  particles. 

The  paper  is  organized  as  follows.  Section  (2)  contains  a  general  discus¬ 
sion  of  the  turbulent  trapping  model,  and  Section  (3)  a  brief  review  of  previous 
numerical  work.  Sections  (4-5)  are  analytic  in  nature  and  provide  the  basic 
formulation  of  the  problem.  In  sections  (6-8)  we  continue  the  analytic  work, 
and  review  the  quantitative  predictions  of  the  turbulent  trapping  model  in  some 
detail.  It  is  in  sections  (9-12)  that  we  embark  on  the  numerical  simulations, 
and  present  our  results  for  comparison  with  the  analytic  work.  The  discussion 
of  these  results  is  centered  on  the  main  predictions  of  the  turbulent-trapping 
model,  concerning  the  structure  of  the  correlation  functions,  and  the  enhance¬ 
ments  of  the  growth  rates. 

We  shall  summarize  our  main  conclusions.  (1)  As  regards  the  correlation 
functions  and  the  consequent  structure  of  turbulent  phase  space,  the  turbulent 
trapping  model  provides  a  surprisingly  accurate  prediction,  especially  in  view 
of  its  semi-qualitative  nature.  (2)  While  there  is  evidence  for  an  enhancement 
of  the  growth  rates  above  the  quasilinear  values,  we  have  not  obtained  a  quan¬ 
titative  verification  of  the  predictions  of  the  turbulent  trapping  model.  In  the 
regime  where  the  model  is  expected  to  be  strictly  valid,  the  enhancements  ob¬ 
served  are  considerably  smaller  than  predicted;  and  when  the  enhancements  are 
closer  to  what  would  be  expected  from  the  turbulent  trapping  model,  they  are 
observed  in  a  regime  where  the  model  is  no  longer  strictly  valid.  We  believe 
however  that  these  results  for  the  growth  rates  are  not  final,  and  this  because 
our  simulations  operated  under  numerical  limitations  which  might  be  removed 
in  future  work. 


II.  GENERAL  DISCUSSION  OF  THE  THEORY. 

The  non- validity  of  the  quasilinear  theory  in  a  regime  where  the  latter  was 
thought  to  be  correct  can  be  understood  on  a  formal  basis,  starting  from  the 
fact  that  the  Vlasov-Poisson  equations  form  a  quadratically  nonlinear  set  for  the 
distribution  function.  Thus  the  standard  techniques,  applicable  to  any  quadrat¬ 
ically  nonlinear  equations,  can  be  applied  to  the  Vlasov-Poisson  equations  as 
well.  These  techniques  are  mainly,  the  BBGKY  hierarchy,  the  iterative  meth¬ 
ods  used  in  the  Soviet  literature  in  order  to  derive  the  Wave  Kinetic  Equation 
(with  the  Random  Phase  Approximation),  and  lastly,  the  Direct  Interaction 
Approximation  (DIA).  The  DIA  equations  were  shown  by  D.  Dubois  and  M. 
Espedal(7),  and  by  others(8),  to  add  new  self-consistency  terms  (or  polariza¬ 
tion  terms)  not  present  in  other  renormalized  turbulence  theories,  theories  that 
were  based  on  the  simplifying  assumption  of  quasi-Gaussian  electric  field  fluc¬ 
tuations.  In  fact,  such  self-consistency  terms  are  easily  found  to  be  inherent 
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in  the  BBGKY  hierarchy  and  in  the  iterative  methods  as  well.  Moreover,  the 
mode-coupling  terms  identified  by  ALP  to  give  a  contribution  of  the  same  order 
as  quasilinear  may  be  shown  to  belong  to  this  class  of  self-consistency  terms. 
The  non-validity  of  the  quasilinear  theory  in  the  regime  i/*  >  7*  is  a  direct 
consequence  of  the  self-consistency  of  the  electric  field:  namely,  the  statistical 
properties  of  the  electric  field  are  self-consist ently  related  to  the  particle  mo¬ 
tions  through  the  Poisson  equation,  this  self-consistency  makes  the  electric  field 
deviate  from  Gaussian  statistics,  and  this  in  turn  leads  to  a  modification  of  the 
diffusion  coefficient  in  the  regime  Vk  >  7*.  This  result  has  to  be  contrasted 
with  the  so-called  stochastic  acceleration  problem  in  which  the  electric  field  is 
externally  given.  In  this  case  the  statistical  properties  of  the  electric  field  may 
be  assumed  Gaussian,  and  the  quasilinear  predictions  follow. 

A  ”  bootstrap”  argument  can  help  to  provide  a  qualitative  description  of 
these  effects  Let  us  assume  that  we  start  off  with  a  rather  particular  organiza¬ 
tion  of  the  turbulent  fields,  one  that  partitions  them  into  a  set  of  wave  packets 
of  different  phase  velocities  and  finite  amplitudes,  which,  for  a  given  phase  ve¬ 
locity,  have  small  overlap  in  real  space,  and,  when  considered  in  velocity  space, 
have  barely  overlapping  trapping  widths  In  this  picture,  the  envelopes  of  the 
wave  packets  are  almost  stationary,  as  the  group  velocity  of  the  electron-plasma 
waves  is  very  small.  The  beam  particles  stream  through  these  almost  immobile 
packets,  to  be  momentarily  trapped  and  scattered  on  their  way. 

If  one  attempts  to  partition  a  turbulent  electric  field  of  given  total  en¬ 
ergy  and  given  total  spectral  width  into  such  a  configuration  of  wave  pack¬ 
ets,  one  finds  that  the  characteristic  spectral  width  of  each  packet  is  of  order 
6k  «  k(vk/wp),  where  Vk  =  ( k2Dql )1^3  is  the  ” classical”  resonance  broaden¬ 
ing  £requency(9)  and  Dql  is  the  quasilinear  diffusion  coefficient,  proportional  to 
\Ek\ 2.  Each  wave  packet  has  a  trapping  width  of  order  (D^/^)1^3?  has  a  length 
in  real  space  of  order  1  /6k,  and  is  such  that  a  particle  will  reside  in  it  for  just 
about  one  trapping  time,  i/jj”1.  One  can  show  that  such  a  configuration  leads 
to  non- Gaussian  statistics  of  the  electric  fields , 

The  fact  that  the  separatrices  of  the  wave  packets  just  barely  overlap  is 
essential  to  our  argument.  Under  such  a  circumstance,  though  in  most  of  the 
phase  space  near  the  wave  packets  particle  motion  is  stochastic,  diffusion  in  a 
regime  so  close  to  the  stochasticity  threshold  need  not  be  quasilinear.  Thus, 
the  local  diffusion  coefficient  may  differ  from  the  quasilinear  one  by  a  factor  of 
order  1(10).  If  one  then  links,  through  Poisson’s  equation,  the  energy  loss  of  the 
particles  to  the  equal  energy  gain  of  the  wave  packets,  one  must  conclude  that 
at  least  initially,  starting  from  this  particular  configuration  of  wave  packets, 
the  electric  field  energy  will  grow  at  a  rate  7*  different  from  the  quasilinear  In 
other  words,  the  motion  of  an  individual  particle  is  sufficiently  stochastic  to  be 
diffusive  over  long  times,  yet  is  sufficiently  coherent  over  a  single  trapping  time 
for  a  modification  of  the  particle’s  wave  emission  rate  to  occur. 

Now,  two  cases  can  obtain.  If  the  growth  rate  is  larger  than  the  trapping 
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frequency,  7*  >  I/*,  the  wave  packet  structure  will  be  destroyed  by  overlap  of 
neighbors  before  the  modified  diffusion  and  wave  emission  can  occur.  In  this 
case,  we  can  conclude  that  the  initial  configuration  of  "barely  overlapping" 
wave  packets  is  irrelevant,  because  it  subsists  for  a  time  much  shorter  than  the 
diffusive  time-step.  The  second  case,  7*  <  i/*,  is  the  one  of  interest  to  us.  In 
this  limit,  particles  are  trapped  and  detrapped  many  times  in  one  e-folding, 
and  the  wave  packets  can  be  thought  of  feeding  quasistatically  on  the  nonlinear 
dynamics  of  the  particles,  at  a  rate  different  from  quasilinear. 

The  "bootstrap”  argument  is  valid  only  if  we  can  answer  two  questions, 
which  apply  in  the  regime  7*  <  i/*:  (1)  how  does  the  wave-packet  configuration 
we  have  chosen  to  discuss  arise  in  the  first  place,  and  (2)  how  does  it  sub¬ 
sist  through  many  e-foldings,  when  individual  wave  packets  must  necessarily 
disappear,  to  be  replaced  by  others? 

To  answer  question  (1),  we  can  argue  that  if  the  modified  growth  rate 
regime  results  in  an  enhancement  of  the  growth  rate,  7 *  >  7^,  then  the  wave 
packet  structure  will  eventually  dominate  any  other.  This  argument,  which 
is  valid  in  the  case  of  a  linear  instability  dominating  all  other  linear  waves, 
is  harder  to  justify  in  this  case  because  the  phenomenon  is  nonlinear  on  a 
microscopic  scale.  We  shall  however  keep  to  this  qualitative  answer. 

Question  (2)  is  intimately  linked  to  the  self-consistency  of  the  electric  field. 
We  have  no  definite  qualitative  answer,  but  can  try  to  rephrase  the  question  in 
terms  of  a  local  energy  exchange  mechanism.  We  note  that  as  the  wave  energy 
grows,  the  trapping  width  of  the  average  wave  packet  increases.  For  a  constant 
spectral  width,  this  means  that  at  a  given  point  in  space  the  number  of  wave 
packets  will  diminish  with  time.  This  suggests  a  local  mechanism,  by  which  one 
wave  packet  can  dominate  over  its  neighbors,  by  diverting  the  particles  on  whose 
energy  they  would  have  grown,  to  the  extent  of  completely  depleting  some  of 
these  neighbors  and  extending  its  trapping  width  to  their  region  of  phase  space. 
Assuming  this  mechanism  exists,  we  have  a  picture  of  a  configuration  slowly 
transforming  itself  over  a  time  scale  comparable  to  the  growth  rate,  through  a 
process  of  competition  between  neighboring  wave  packets.  Now,  such  a  process 
is  a  coupling  of  modes,  between  waves  differing  by  6k  ~  kiy^/wp)  and  thus, 
the  presence  of  mode  coupling  terms  in  the  analytic  theory,  imposed  by  self- 
consistency,  is  consistent  with  this  picture  of  a  self-sustaining  turbulence  of 
wave  packets. 

Our  picture  is  analogous  to  the  self-sustaining  clump  turbulence  described 
by  Dupree(u),  but  with  the  important  difference  that  in  our  case  the  plasma 
waves  display  a  well-defined  dispersion  relation  w  =  imposed  by  the  linear 
bulk  plasma,  and  which  is  absent  in  the  formulation  of  clump  turbulence.  We 
shall  see  that  this  results  in  correlation  functions  of  the  distribution  function 
which  extend  over  many  wave  periods,  in  sharp  contrast  to  the  clump  correla¬ 
tion  function,  whose  coherence  length  is  limited  at  most  to  one  wave  period. 
Notwithstanding  these  differences,  and  extending  the  notion  of  "clump”  to  the 
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correlation  function  in  the  beam-plasma  instability,  one  might  reformulate  the 
wave-packet  picture  of  the  preceding  paragraphs  by  saying  that  the  modifica¬ 
tion  of  the  growth  rate,  in  the  clump  picture,  results  from  enhanced  Cerenkov 
emission  by  the  clumps,  which  radiate  into  a  background  of  adiabatic  beam 
particles  (see  also  Ref.(,12,13)). 


III.  PREVIOUS  NUMERICAL  WORK. 

While  simulations  of  the  beam-plasma  instability  are  nothing  new  (,14,15), 
performing  such  simulations  under  conditions  where  quasilinear  theory  can  be 
verified  without  ambiguity  is  demanding,  both  in  computational  resources  and 
in  the  care  with  which  the  parameters  of  the  simulation  must  be  chosen.  This  is 
because  behind  its  apparent  simplicity,  the  weak  beam-plasma  instability  hides 
several  time  scales,  which  must  be  correctly  resolved  for  quasilinear  theory  to 
be  applicable  (,16,17)  The  same  applies  to  the  turbulent  trapping  model  which 
assumes  constraints  similar  to  those  of  quasilinear  theory. 

To  our  knowledge,  the  first  numerical  experiments  testing  quasilinear  the¬ 
ory  in  its  full  regime  of  validity  are  due  to  J,C.  Adam(18).  These  simula¬ 
tions  exhibited  the  formation  of  discontinuous  electric  field  spectra,  ascribed 
to  mode-coupling  effects,  and  showed  an  enhancement  of  growth  rates  over  the 
quasilinear  values,  but  unfortunately  were  not  sufficiently  detailed  for  quanti¬ 
tative  conclusions  about  the  turbulent  trapping  regime  to  be  reached.  A  set  of 
older  and  somewhat  complementary  numerical  experiments  are  those  due  to  A. 
Bakai  and  Yu.  Sigov  (19),  who  emphasized  in  their  study  the  qualitative  obser¬ 
vation  of  phase  space  structure,  without  however  providing  more  quantitative 
descriptions  of  the  phenomena. 

In  the  present  work,  we  have  tried  to  proceed  well  beyond  the  results  of 
these  previous  simulations,  by  a  more  systematic  measurement  of  growth  rates, 
by  the  quantitative  description  of  phase  space  through  the  use  of  correlation 
functions,  and  finally,  by  the  direct  comparison  of  numerical  results  with  the 
analytic  predictions  of  the  turbulent  trapping  model. 


IV.  PHYSICAL  PROBLEM  AND  APPROXIMATIONS. 

The  plasma  is  one-dimensional  and  consists  of  one  species,  electrons  moving 
against  an  immobile  neutralizing  background  of  ions.  A  sketch  of  the  initial, 
unperturbed  distribution  function  is  shown  in  Fig.(l). 

In  our  simplified  model,  the  bulk  plasma  is  completely  cold  and  is  rep¬ 
resented  by  linear  fluid  equations  for  the  perturbed  bulk  velocity  and  density. 
The  beam  plasma  on  the  other  hand,  is  described  by  the  exact  Vlasov  equation, 
with  no  other  approximation  in  the  numerical  simulations,  than  those  imposed 
by  the  finite  differencing  method  we  have  chosen.  This  numerical  approach, 
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in  which  the  bulk  is  ’’streamlined”  so  as  to  simply  support  linear  plasma  os¬ 
cillations,  is  valid  for  problems  where  the  detailed  beam  dynamics  are  most 
important.  Modelling  the  bulk  as  cold  and  linear  then  results  in  a  considerable 
or  even  vital  economy  of  mesh  size. 

Let  the  total,  averaged  spatial  density  of  the  electrons  be  no,  with  this 
density  split  between  the  bulk,  nop,  and  the  beam,  non,  nop  4*  Kqb  =  fto*  The 
beam  is  a  small  perturbation  to  the  bulk,  n op  nop.  We  model  the  bulk 
evolution  by  the  coupled  continuity  and  momentum  equations: 


—Tip  =  -Tlop—Up,  (1) 

%-uP  =  ±E(x,t)  (2) 

eft  m 

where  hp  and  up  are  the  linearized  fluid  variables.  The  beam  is  described 
by  a  distribution  function,  /b(x,v,£)  which  evolves  according  to  the  Vlasov 
equation: 


d  M  a  £  q  d  £ 

■kiJb  +  v— /b  H - q~/b  =  0, 

at  ox  mov 


(3) 


with  the  normalization  of  the  unperturbed  distribution  function: 


dv  fB(v) 


=  no  B, 


(4) 


The  electric  field  E(x ,  t)  is  determined  self-consist ently  from  Poisson’s  equation: 


E(x,t )  —  ^nop  +  nP(x,t)  +  J dvfB(x,v,t)  -  noj  ,  (5) 

where  we  have  subtracted  from  the  electron-charge  source  terms  the  neutralizing 
charge  of  the  immobile  background  ions,  qtiQ. 

We  proceed  to  normalize  eqs.(l-  5).  We  already  have  a  characteristic  time- 
scale,  determined  by  the  plasma  frequency: 


0 fp 


me o  ’ 


(6) 


but,  because  the  bulk  is  cold,  we  cannot  choose  its  thermal  velocity  as  a  refer¬ 
ence  velocity,  nor  the  Debye  length  A]?  =  v^/wp  as  a  reference  length,  as  both 
are  zero.  Rather,  we  introduce  a  reference  velocity  vr  which  is  of  the  order  of 
the  beam  velocities  (Fig.(l)).  The  reference  length  is  then  A r  =  vr/ujp ,  and 
we  define  the  normalized  variables  as: 
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t'  =  uPt,  x'  =  x/XR ,  v'  =  v/vR, 


(7) 


E'  =  — ?-■  , 

T7ia;pVR 


up  =  up/vr,  n'p  =  np/n0p, 


(8) 


/'  = 


— /b, 

«ob 


Dropping  the  primes,  eqs.(l-  5)  then  become: 


(9) 


a  a 

atnp  “  -&*'• 

I  *'■  =  E- 


(10) 

(ii) 


(12) 


— =  iip  np  -f  jRp  /  dt;  —  1^  .  (13) 

with  the  normalization: 

J  dx  J  dvf(x,  vy  t)  =  L  (14) 

where  L  is  the  system  length,  and  where  we  have  defined  the  parameters, 


d  _  no* 
—  ? 
no 


Rp  =  ^2L  =  l-RB, 
no 


(15) 


which  gauge  the  relative  bulk  and  beam  densities,  with  iip  1.  IfJRp—0  ex¬ 
actly,  Eqs.(10,ll,13)  predict  linear  plasma  oscillations,  with  a  plasma  frequency 
ojp  =  1.  For  convenience,  we  shall  drop  all  primes  in  referring  to  the  normalized 
equations. 

For  eqs.(10-  13),  we  can  define  the  normalized  energies  of  the  electric  field, 
of  the  bulk  plasma  and  of  the  beam  plasma  particles.  These  are,  respectively: 


to£  =  i  J  E2(x,t)  dx, 

(16) 

wP  =  i RP  J u2p(x,t)  dx, 

(17) 

wB  =  | Rb  J  J  v2f(x,v,t)  dx  dv. 

(18) 
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and  with  these  definitions,  conservation  of  the  total  energy  of  the  system  follows: 

liWtot  =  ^(we  +  u>p  +  u>b)  =  0,  (19) 


V.  LINEAR  GROWTH  RATES  AND  QUASILINEAR  THEORY. 
For  a  very  weak  beam,  Rb  1,  eqs.(10-  13)  have  the  dispersion  relation: 

u>  =  ±u>k ,  Ljk  =  sgn(fc),  |u>fc|  =  1,  (20) 

so  that  the  phase  velocities  are  t/*  =  Uk/k  =  ±l/|fe|.  With  this  particular 
definition,  the  waves  on  the  branch  all  propagate  to  the  right,  and  those 
on  the  branch  all  to  the  left.  The  linear  growth  rate  of  the  waves  is  given 
by: 


7*  =  <  f'(vk)  >„  (21) 

where  <  f(v)  >,  is  the  space-averaged  distribution  function  of  the  beam  (see 
Appendix  A  for  a  discussion  of  ensemble  and  spatial  averages).  For  the  finite- 
length  system  of  the  numerical  simulations,  the  mode  structure  is  discrete, 
with  wave  numbers  k  =  2nn/L ,  n  =  1,2,. ...  The  electric  field  has  the  modal 
expansion: 


E(x,*)  =  ^£+(t)ei(fc— (22) 

k  k 

where  we  made  explicit  right  and  left-propagating  waves.  With  the  beam  dis¬ 
tribution  function  confined  to  positive  velocities,  only  the  right-propagating 
modes  E£  (t)  are  resonantly  driven  by  the  beam  particles,  and  their  amplitudes 
eventually  dominate  over  the  left-going  waves  (t). 

The  correlation  time  of  the  electric  fields,  defined  as  the  width  of  the  two- 
point  correlation  function,  is  given  roughly  by: 

_  2?r _ 2tt 

Tc~A (u/fc-fcv)  ~  t;  Afc*’  (  } 

where  Aks  refers  to  the  total  width  of  the  turbulent  spectrum  excited  by  the 
beam.  In  the  usual  formulation  of  quasilinear  theory,  it  is  understood  that  the 
correlation  time  is  very  short,  and  that  we  have  the  ordering: 


re  <  7fc  1  »  t</>> 


(24) 
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where  r</>  is  the  characteristic  time  for  the  evolution  of  <  f(v)  >,.  An 
equivalent  way  of  stating  Eq.(24)  is  through  the  so-called  O’Neil-Malmberg 
parameter  77  (20),  which  can  be  written,  with  our  choice  of  normalizations: 

77  =  Rb{vb/&vb)3,  (25) 

In  this  expression,  vb  is  the  mean  phase  velocity  of  the  unstable  waves.  With 
this  definition,  the  condition  77  -C  1  can  be  shown  to  be  equivalent  to  Eq.(24). 

If  we  use  the  standard  formula  for  the  quasilinear  diffusion  coefficient  in 
the  resonant  velocity  region  (see  for  instance(21)),  we  find: 

Dqt(v)  =  (26) 

where  =  1/v  and  where  L  is  the  total  length  of  the  system.  The  overbar 
average  is  defined  in  Eq.(51)  of  Appendix  A:  it  is  meant  to  reconcile  the  star 
tistical  properties  of  a  single,  long  system  with  those  of  an  ensemble  of  such 
systems, 


VI.  THE  EQUATION  FOR  THE  CORRELATION  FUNCTION. 

The  starting  point  for  the  turbulent  trapping  model  is  the  2-point,  1-time 
correlation  function,  which  is  a  measure  of  the  phase-space  ”  graininess”  of  the 
beam  distribution  function.  It  is  defined  by: 


C(x_,t/_,v+,t)  =  <  6f(zuvut)  6f(x2,V2,t)  >,  (27) 

where  x_  =  Xi  —  x2,  =  Vi  —  t>2,  v+  =  (vi  4*  v2)/2 ,  and  where  6f  is  the 

fluctuation  of  the  distribution  function  about  the  averaged  distribution  6f  = 
/—  <  /  >.  Implicit  in  Eq.(27)  is  the  assumption  that  the  turbulence  is  spatially 
homogeneous,  so  that  C  does  not  depend  on  x+  =  (xi  +  x^j/2.  Furthermore, 
C  is  a  slowly-varying  function  of  t/+  (with  variation  on  the  scale  of  the  total 
beam  width),  so  that  we  shall  keep  this  variable  a  tacit  parameter  by  writing 
C  =  C(x_, u_,  t)  in  most  places. 

On  the  basis  of  a  number  of  simplifying  assumptions  done  in  the  spirit  of 
the  turbulent  trapping  model,  one  can  derive  (6)  a  Fokker-Planck  equation  for 
C(x_,v_,t)  which  has  the  form: 

(dt  +  v-dx_  -  2  (1  -  cosk+x-)D(v+, t)d%_^  C(x-,v~,t)  = 

2  cosk+x-  D(v+,t)  ( dv+  <  f  >)2, 


(28) 
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where  =  l/v_|_  is  the  wavenumber  resonant  with  the  average  velocity  v+.  An 
important  feature  of  Eq.(28)  is  that  the  diffusion  term  £>(u+,£)  is  left  undeter¬ 
mined.  It  is  defined  as  the  one-particle  diffusion  coefficient,  for  an  ensemble  of 
particles  starting  at  t/(£)  =  and  executing  a  diffusive  motion  in  the  turbulent 
fields  until  a  later  time  tf  =  t  -f  T.  That  is: 

D(v+<>t)  =<  ^  2 ^ 

where  Av  —  v(t  +  T)  —  v(t).  Provided  the  time  scales  are  well  separated,  there 
will  exist  a  diffusive  regime  rc  <T  7^X>  such  that  Eq.(29)  provides  a  value 
of  D  independent  of  T.  Now,  if  quasilinear  theory  is  valid,  then  from  Eq.(29)  we 
must  necessarily  obtain  D(v+,£)  =  Dql(v+it),  where  Dql  is  given  by  Eq.(26). 
However,  Eq.(28)  allows  for  a  more  general  type  of  diffusive  motion.  In  fact,  the 
only  assumption  we  retain  is  that  D(t/+,f)  grows  as  the  square  of  the  resonant 
wave  amplitude,  so  that  £)(t;+,t)  «  exp(27*+t).  Because  7^  ^  r</>>  this 
means  that  the  entire  right  hand  side  of  Eq.(28)  has  this  exponential  dependence 
as  well,  to  within  the  slow  modulation  by  (5V+  <  /  >)2,  which  occurs  on  the 
T<f>  time  scale.  Thus  the  time  derivative  on  the  left  hand  side  of  Eq.(28)  is  of 
order  dt  «  27*+ . 

Let  us  ”  freeze”  the  value  of  D  at  a  given  time,  and  define  an  instantaneous 


trapping  frequency  and  an  instantaneous  trapping  velocity: 

% 

Vtt  =  (k+D(v+))113  ,  A vtt  =  {D(v+)/k+)1/3 ,  (30) 

If  we  then  normalize  variables  in  terms  of  these  trapping  parameters: 

u  =  v-/2l/zAvtt,  £  =  fc+x_,  r  =  21^3i/«  t  (31) 

C(x. ,  v.)  =  22'3 Av?t  ( 8V+  <  f  >)2  H(t,  u)  (32) 

we  obtain  the  normalized  equation: 

(  6r  +  -  (1  -  cosf)^  )  H(Z,  u)  =  cost,  (33) 


where  6r  =  0(jk+/^u)  symbolizes  the  importance  of  the  time  derivative:  Eq.(33) 
is  in  fact  rigorous  only  provided  the  time-derivative  term  6r  can  be  neglected, 
which  requires  7fc+  i /«. 

Eqs.(28)  or  (33)  incorporate  two  basic  features  of  the  turbulent  motion 
of  the  beam  particles.  The  first  is  that  as  x_  — ►  0  and  v_  — ►  0,  C(x.,u_) 
becomes  large,  reflecting  the  fact  that  particles  initially  very  close  in  phase  space 
remain  correlated  for  a  long  time,  and  do  not  random-walk  independently  in 
velocity  space(u).  The  second  feature  of  Eq.(28)  is  that  the  diffusion  coefficients 
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are  periodic  in  i.,  with  a  period  equal  to  that  of  the  wave  with  which  the 
mean  velocity  of  the  particles  is  resonant.  Again,  this  is  linked  to  the  fact 
that  particles  initially  an  integral  number  of  wavelengths  apart  will  undergo 
strongly  correlated  motion  This  second  property  is  only  approximate,  it  will 
be  discussed  m  more  detail  m  the  next  section. 


VII.  SOLUTION  IN  THE  TURBULENT  TRAPPING  REGIME. 

The  regime  of  turbulence  is  determined  by  the  relative  importance,  in 
Eq.(33),  of  the  time-derivative  term  ST  =  0(7*+/l/tt)*  This  distinguishes  two 
regimes.  In  the  first,  the  growth-rate  dominates  the  trapping  frequencies, 
jql  >>  vti ,  and  in  this  case  we  can  neglect  the  d^_  term  in  Eq.(28).  In  effect, 
the  fast  growth  rate  washes  out  any  harmonic  strucuture  linked  to  the  electric 
field  turbulence  (apart  from  the  overall  amplitude).  G.  Laval  and  D  Pesme(6) 
have  shown  that  in  this  regime  one  recovers  the  quasilinear  growth-rate. 

The  ” turbulent  trapping”  regime  exists  in  the  opposite  limit: 

79t  «V't,  (34) 

In  this  limit,  we  neglect  6r  in  eq.(33),  and  we  can  say  that  the  modes  grow 
”quasistatically”  when  their  slow  growth  is  measured  on  the  now  relatively 
short  time-scale  of  the  turbulent  trapping,  i^J1.  With  6T  ignorable,  Eq.(33)  is 
rigorous  and  has  an  exact  series  solution(6),  given  by: 


H(£,u)  =  ^2  Hn{u)  exp(tn(£  -  sin£)),  (35) 

n— — oo 

where  the  sum  excludes  the  term  n  =  0,  and  with: 

ff»(«)  =  p(u  N 1/3)  (36) 

with  N  =  |n|,  crn  =  sgn(n),  and  with  the  Airy-related  function: 

P(z)  =  -  [°°  exp  (— T3/ 3  -  izT )  dT  (37) 

*  Jo 

Graphs  of  if(^,u),  obtained  from  the  analytic  series  (35),  are  shown  in 
figs.  (2).  A  salient  property  of  this  approximate  solution  is  that  H(£,u)  has  a 
logarithmic  divergence  as  £,  u  — ►  0,  with  the  asymptotic  expression: 


1  3 

W  ^73  1o6  ^  _  22/3^u  +  21/3U2)  ’ 


(38) 
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This  divergence  reflects  the  property  required  of  Eq.(28),  that  points  of  phase 
space  initially  very  close  remain  correlated  for  long  times.  A  similar  struc¬ 
ture  arises  in  the  solution  of  T.H.Dupree’s  equations  for  turbulence  in  a  bulk 
plasma(11). 

The  series  solution,  Eq.(35),  diverges  as  f ,  u  — ►  0  because  we  have  neglected 
the  <5r  term  in  Eq.(33),  a  term  reflecting  the  importance  of  the  time  derivative  in 
Eq.(28)  for  small  separations  and  t;_.  Indeed,  we  expect  (35)  to  break  down 
when  the  terms  ud^  and  (1  —  cosf)c^  in  Eq.(33)  are  comparable  to  ST.  If  we 
assume  that  «  du  «  1,  in  the  normalized,  approximate  solution  if(£,u),  then 
the  <5r  term  becomes  important  when  u  <  8r  and  £  <  8r  -  In  physical  units, 
the  conditions  for  the  validity  of  (35)  are  then  |v_|  7*+//c+  and  |fc+x_| 

(7fc+Mt)1/2. 

While  an  analytic  solution  with  dt  ^  0  has  not  been  found,  we  can  obtain 
the  special  value  of  C(0, 0)  very  simply,  by  letting  x_,  v_  — ►  0  in  Eq.(28),  which 
eliminates  all  but  the  first  term  on  the  left.  We  have,  with  dt  =  27*+ : 

C( o,  0)  =<  (« 5/)2  >=  {dv+  <f>)2,  (39) 

7*+ 

an  expression  which  assumes  that  7^ 

We  should  also  stress  that  the  strict  periodicity  of  C(x_,  v_)  as  a  function 
of  x_  is  an  approximation  which  results  from  assuming  that  particles  at  t;  =  t/+ 
feel  the  single  wave  number  =  &+  =  l/v+.  To  be  consistent  with  the  spirit 
of  the  turbulent  trapping  model,  we  must  in  fact  assume  that  the  particles  see 
an  entire  wave  packet,  still  centered  about  k+y  but  of  width  6ktt  =  The 

maximum  correlation  length,  will  be  of  the  order  of  the  length  of  this  wave 
packet: 

Itt  =  1  /k+vtt  =  l/k^D1'3,  (40) 

Qualitatively,  we  expect  the  series  solution  (35)  to  be  modified  by  an  enve¬ 
lope  function  of  width  lti ,  such  that  successive  peaks  of  C(x-,v~)  diminish  in 
amplitude,  tending  to  zero  for  |x_|  ^  Itt-  As  we  shall  see,  this  behavior  is 
numerically  confirmed. 

By  inspection  of  Fig.  (2a),  we  can  see  that  the  total  width  of  the  profile  of 
H(u,  0)  is  about  3.  This  width  in  t/_  is  then  approximately  3  X  21/3  Avtt  ~ 
AAvtt.  Now,  for  Eq.(28)  to  be  valid,  the  width  in  t/_  of  the  correlation  function 
must  be  small  compared  to  the  beam  width.  This  imposes  the  reasonable 
trapping  width  constraint: 


4Avtt  <C  A vbj 


(41) 
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where  Avp  is  the  beam  width  over  which  waves  axe  excited.  In  other  words,  any 
turbulent  structure  must  remain  small,  and  it  is  in  this  sense  that  the  turbulent 
trapping  model  remains  within  weak  turbulence  theory. 


VIII.  ANALYTIC  GROWTH  RATES  IN  THE  TURBULENT  TRAP¬ 
PING  REGIME. 

By  using  Poisson’s  equation,  we  can  express  the  growth  rate  7*  in  terms  of 
a  velocity  integral  of  the  correlation  function.  G  Laval  and  D.  Pesme(6)  have 
used  this  expression  to  find  the  growth  rate  in  the  turbulent  trapping  regime, 
by  what  might  be  termed  a  ” bootstrap”  method.  First,  C(z_,  v_)  is  derived 
from  the  series  solution  in  Eq.(35).  With  Poisson’s  equation,  this  results  in  an 
expression  where  7*  is  an  explicit  function  of  Z>,  the  undetermined  diffusion 
coefficient.  It  is  then  found  that  by  choosing  a  particular  form  for  7*  and  D, 
global  energy  conservation  is  satisfied  while  the  functional  relation  between  7* 
and  D  remains  true.  The  expressions  for  7*  and  D  are  simply: 


7*  =  2.2  72',  (42) 

D{ v)  =  2.2  Dql{v),  (43) 

where  72'  and  Dql(y)  are  given  by  eqs.(21)  and  (26).  The  numerical  factor  of 
2.2  comes  from  the  summing  of  a  series  of  harmonics  in  Eq.(35),  in  the  form 
of  a  series  of  Bessel  functions.  The  striking  feature  of  eqs. (42,43)  is  that  the 
enhancement  of  the  growth  rate  and  diffusion  coefficient  is  this  simple  numerical 
factor,  independent  of  k  and  of  the  level  of  turbulence,  provided  vtt  7*. 

Analytic  work(4),  done  in  the  complementary  regime  utt  <  7*,  suggests  a 
refined  threshold  for  turbulent  trapping,  with  the  strong  inequality  replaced  by 
the  condition: 


>  5  7fc, 


(44) 


The  work  of  Ref.(4)  also  suggests  that,  starting  from  initial  conditions,  the 
enhancement  of  the  growth  rate  in  the  turbulent  trapping  regime  will  occur 
only  after  a  substantial  time  lapse,  that  is  after  at  least  a  few  e-foldings  of  the 
electric  field.  This  is  because  for  eqs.  (42,43)  to  be  valid,  the  non-resonant  har¬ 
monics  of  the  Fourier  components  of  the  electric  field  must  be  in  a  quasistatic 
equilibrium  with  the  fundamental,  resonant  harmonics.  At  t  =  0,  the  nonreso- 
nant  components  are  rigorously  zero,  and  thus  a  finite  time  must  elapse,  before 
the  mode-coupling  mechanisms  can  build-up  the  hierarchy  of  modes,  mutually 
scattering  into  each  other,  which  contribute  to  the  enhancement  of  the  growth 
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of  the  fundamental.  If  N  is  the  number  of  e-foldings  undergone  by  the  electric 
field,  we  shall  require: 


J\T«  3  -5, 

for  a  regime  of  enhanced  growth  rates  to  manifest  itself. 


(45) 


DC.  THE  NUMERICAL  SCHEME. 

The  Vlasov  solver  for  Eq.(10)  uses  a  finite  difference  scheme  developped  by 
J.P  Boris  and  D.L.  Book(22)  ,  the  so-called  ” flux-corrected  transport”  method, 
which  exhibits  a  particularly  low  diffusivity,  while  maintaining  good  numerical 
stability  even  with  very  sharp  density  profiles  in  phase  space.  This  scheme  is 
combined  with  a  split-step  method  devised  by  C.Z.  Cheng  and  G.  Knorr(23), 
in  which  transport  of  the  distribution  function  alternatively  proceeds  in  x  and 
in  v . 

The  distribution  function  for  the  beam  particles,  /(x,i/),  is  defined  on  a 
mesh  with  Nx  points  in  the  spatial  dimension  x  and  Nv  points  in  velocity  r, 
with  spacings  Ax  and  Av  respectively.  The  total  system  length  i sL  =  Nx  Ax, 
with  boundaries  in  x  at  Xmin  —  0,  xma*  =  L.  The  boundaries  in  velocity  are  at 
v  —  vmin  and  v  =  vmax.  The  spatial  boundary  conditions  are  periodic,  while  at 
v  =  vmin  ?  vmax  we  impose  /  =  0,  a  good  approximation  provided  /  is  already 
very  small  some  distance  from  the  velocity  boundaries.  We  advance  the  system 
in  time  with  a  time-step  At,  usually  chosen  so  that  up  At  =  0.2. 

The  bulk  equations  are  solved  independently  on  the  one-dimensional  spatial 
mesh,  using  fast  Fourier  transforms  in  a  leapfrog  scheme.  Beam  and  bulk 
equations  are  then  coupled  through  Poisson’s  equation,  also  solved  with  Fourier 
transforms. 

Because  our  simulations  are  aimed  at  exploring  the  initial  turbulence  cre¬ 
ated  in  the  beam-plasma  instability,  we  think  of  the  beam  distribution  function 
as  a  tunable  source  of  free  energy,  and  we  define  its  shape  to  obtain  a  con¬ 
stant  initial  growth  rate(19)  over  a  large  band  of  phase  velodties(figs.(3a,b)). 
This  is  done  by  piecing  together  simple  functions  of  velocity  (see  Appendix 
C).  Our  choice  of  an  initial  distribution  function  is  somewhat  artificial,  but  it 
has  the  great  advantage  of  avoiding  trapping  effects  which  arise  from  the  early 
dominance  of  a  single  mode. 

In  Table  I,  we  list  the  principal  parameters  of  the  two  simulations  to  which 
we  shall  refer  in  the  following  sections.  In  the  table,  we  indicate  how  well  the 
validity  conditions,  summarized  by  eqs. (54-58)  of  appendix  B,  are  satisfied. 

We  initially  perturbed  the  beam  distribution  function  and  the  bulk  plasma 
density  with  a  purely  spatial  modulation,  so  as  to  obtain  an  initial  electric  field: 
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k , 

E(x, 0)  =  ^  efcsin (kx  +  fa),  (46) 

k=ki 

where  the  phases  <f>k  are  chosen  randomly  from  one  mode  to  the  next,  and 
where  the  coefficients  e*  are  defined  so  as  to  provide  an  electric  field  with  a 
given,  smooth  amplitude  spectrum.  The  modes  in  ki  <  k  <  are  restricted 
to  lie  in  the  unstable  part  of  the  spectrum.  (  Fig.  (3a)  ). 

The  initial  standing  wave  of  Eq.  (46)  splits  into  right-propagating  and  left- 
propagating  waves,  with  only  the  former  amplified  by  the  beam-plasma  inter¬ 
action.  Neglecting  ballistic  effects,  from  Eq.(22)  we  have  the  linear  behavior: 

1  h 7 

E(x ,  t)  =  -  ^2  e*  [®n(k*  —  t  +  <}>k)elht  +  sin(fcx  +  t  +  <£*)],  (47) 

2  k=kx 

XL  THE  NUMERICAL  CORRELATION  FUNCTION. 

In  this  section,  we  consider  the  numerical  results  obtained  in  Simulation  1 
(Table  I),  for  the  correlation  functions  of  the  beam  distribution  function.  The 
initial  distribution  function  used  in  the  simulation,  the  resulting  linear  growth 
rates  and  the  phase  velocities  of  the  modes  initially  excited  are  shown  in  figs.  (3). 
The  initial  growth  rate  was  yL  =  2.19  x  10-3  (obtained  with  a  relative  beam 
density  Rb  =  nE  /no  =  2x  10"3),  and  the  modes  were  excited  in  the  range  of 
phase  velocities  1.15  <  Vk  <  2.55  (79  modes  in  0.392  <  k  <  0.871)  For  the  form 
of  the  initial  excitation,  we  chose  a  nonuniform  distribution  of  mode  amplitudes, 
with  the  dependence  on  wave  number  e*  ~  l//:3.  This  choice  insured  that  the 
same  amount  of  resonance  overlap  obtained  for  all  modes  of  the  spectrum. 
The  additional  advantage  of  this  uneven  distribution  of  amplitudes  is  that  it 
permitted  us,  at  any  one  instant  in  time,  to  measure  the  correlation  function  for 
an  entire  range  of  trapping  widths  At by  scanning  v  =  across  the  width 
of  the  beam. 

If  we  consider  a  mode  in  the  center  of  the  spectrum,  say  with  v*  =  1.75  ( 
k  =  0.564  ),  we  find  that  initially  Vu/~fL  «  4.5,  and  at  the  end  of  the  simulation 
vtthb  «  9-  Thus  according  to  Eq.(44)  the  center  of  the  spectrum  was  in  the 
turbulent  trapping  regime  from  the  very  onset  of  the  instability.  However,  Sim¬ 
ulation  1  proceeded  for  only  about  one  e-folding,  to  t  =  500,  an  evolution  time 
which  is  probably  insufficient  for  the  manifestation  of  enhanced  growth  rates 
(Eq.(45)).  The  discussion  of  the  latter  effects  is  deferred  to  the  next  sections, 
where  we  discuss  them  in  relation  to  a  longer  simulation  (3.2  e-foldings). 

The  evolution  of  the  total  field  energy  WE{t)  is  shown  in  Fig. (4).  The 
oscillating  character  of  wE{t)  comes  from  the  interference  of  the  growing  right- 
propagating  waves  with  the  left-propagating  waves  of  constant  amplitude. 


-18- 


The  final  space-averaged  distribution  function, <  f(v)  >My  is  shown  in 
Fig. (5),  where  it  is  compared  on  the  same  graph  to  the  initial  distribution 
function.  It  can  be  seen  from  this  figure  that  little  modification  of  <  f(v)  >s 
has  occured,  indicating  that  at  this  point  t</>  >  7”1.  By  looking  at  a  cross 
section  of  /  at  fixed  x,  as  is  done  in  fig(6),  we  can  also  see  that  the  amplitude 
of  fluctuations  in  the  non-averaged  distribution  function  has  remained  small. 

So  that  we  may  obtain  a  more  general  view  of  the  phase  space  organization 
of  the  beam,  in  Fig. (7)  we  show  the  contour  lines  of  the  final  distribution  func¬ 
tion,  /(x,  v)  at  t  =  500,  for  the  entire  length  of  the  simulation,  0  <  x  <  1024. 
We  shall  now  interpret  this  figure.  We  first  note  that  at  t  =  500,  the  wave 
spectrum  has  remained  almost  entirely  confined  to  the  initial  range  of  phase 
velocities,  1.15  <  v*  <  2.55.  Thus  we  expect  trapping  to  occur  within  this 
velocity  range,  and  the  motion  of  particles  to  be  adiabatic  beyond  it,  with  how¬ 
ever  some  extension  of  the  trapping  region  above  and  below  the  range  of  excited 
phase  velocities,  due  to  the  finite  trapping  widths  of  the  waves.  This  is  indeed 
suggested  by  the  contours  of  Fig. (7).  For  v  >  2.8,  the  distribution  function  is 
only  slightly  perturbed  by  the  wave  turbulence,  while  in  the  resonant  or  near- 
resonant  range  2.8  >  v  >  1.1  the  contours  display  a  "graininess”,  consisting 
of  closed  or  almost  closed  patterns  or  ” grains”,  suggestive  of  the  formation  of 
small  plateaus  in  the  distribution  function.  When  they  are  examined  individ¬ 
ually,  these  patterns  have  dimensions  in  x  and  v  which  qualitatively  scale  as 
expected.  At  large  velocities  (say  v  =  2.3),  where  the  resonant  wave  amplitudes 
are  large  (4Av«  =  0.24),  the  circles  or  semicircles  of  the  contours  are  wide  in 
velocity.  At  low  velocities  (say  v  =  1.3),  the  resonant  wave  amplitudes  are  small 
(4Avtt  =  0.15)  and  so  is  the  size  of  the  "grains”.  The  spatial  arrangement  of 
the  "grains”  also  reflects  the  properties  of  the  resonant  wave:  one  can  verify 
that  they  have  a  spacing  roughly  equal  to  the  wavelength  of  the  wave,  which  is 
A+  =  27r /fc+  =  27 rv+. 

To  further  clarify  these  observations,  we  show  in  Fig.  (8)  a  perspective  plot 
of  a  small  portion  of  the  distribution  function,  extracted  from  the  region  which 
is  labelled  "Area  1”  in  Fig. (7).  This  picture  displays  the  plateaus,  which  were 
suggested  by  the  contour  plot  in  Fig.  (7),  in  a  more  convincing  manner,  and  it 
is  from  these  plateaus  that  we  infer  the  existence  of  locally  trapped  particles. 
A  striking  feature  of  the  figure  is  that  phase  space  appears  very  regular  when 
it  is  seen  on  this  small  scale,  a  central  prediction  of  Eq.(28).  As  noted  above  in 
Section  5,  this  regularity  is  only  approximate,  and  should  break  down  beyond 
the  length  of  a  resonant  wave  packet,  a  length  which  we  estimated  to  be  of  order 
ltt  =  \/k+vtt(  Eq.(40)).  For  the  "reference”  velocity  v+  =  1.75,  this  length  is 
ltt  ^  90  at  the  end  of  the  simulation.  In  Fig.  (9)  we  blow  up  the  longer  region 
which  is  labelled  "Area  2”  in  Fig. (7),  and  display  it  in  perspective.  This  picture 
extends  over  a  length  in  x  of  160  units,  almost  three  times  the  length  of  Area  1. 
The  figure  shows  that  when  seen  over  this  larger  space  scale,  the  alignement  of 
successive  plateaus  loses  its  regularity,  indeed  over  a  length  which  is  comparable 
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to  ltt  «  90. 

The  discussion  given  above  is  qualitative  and  visual  in  nature.  To  probe  the 
structure  of  Fig.  (7)  in  a  quantitative  manner,  we  need  a  more  exact  numerical 
tool,  and  such  a  tool  is  precisely  the  correlation  function  of  Eq.(27).  In  the 
simulation,  we  evaluated  C(x_,  v_)  by  performing  a  spatial  average  over  the 
length  of  the  system.  Thus  for  given  values  of  x_  =  t_Ax,  v_  =  j-Av  and 
v+  =  Vjnin  +  (jf+  —  l)Av,  we  calculated  the  sum; 

1  Nm 

C(x_,v_,v+)  =  - ■  fi+i.j++j (48) 

iV*  i= l 

where  fij  is  the  distribution  function  at  the  mesh  point  (i,  j). 

To  make  the  link  with  the  analytic  solution  given  by  Eq.(35),  we  needed  to 
know  the  value  of  the  diffusion  coefficient  in  eqs.(30).  The  simplest  approach 
was  to  assume  that  D  =  Dq*  and  obtain  Dql  from  Eq.(26).  In  so  doing,  we 
ignore  any  possible  increase  in  diffusion,  as  would  be  indicated  by  Eq.(42).  This 
is  justified  to  the  extent  that,  as  mentioned  above,  we  do  not  expect  the  system 
to  have  had  the  time  develop  enhanced  diffusivity  within  a  single  e-folding. 
Furthermore,  should  D  be  enhanced  by  the  amount  suggested  by  Eq.(42),  this 
enhancement  will  lead  to  a  quite  modest  broadening  of  the  trapping  width,  with 
A Vtt  =  ( D/k )x/3  =  1.3Av*t.  From  the  discussion  that  follows,  it  will  be  seen 
that  the  observation  of  such  a  moderate  broadening  is  difficult  to  obtain  from 
the  inspection  of  the  numerical  correlation  function.  Thus,  a  less  ambiguous 
detection  of  enhancement  is  to  be  obtained  from  the  direct  measurements  of 
the  growth  rates,  to  be  presented  in  the  next  section. 

We  first  consider  the  structure  of  the  correlation  function  in  velocity  space, 
by  plotting  C(0,  v_),  for  a  fixed  value  of  the  mean  velocity  v+.  Thus  in  Fig. (10), 
we  compare  the  numerical  results  for  C  (0,  v_ )  with  the  predictions  of  the  series 
solution,  Eq.(35),  for  the  mean  velocity  v+  =  1.745(Av«  =  0.034).  We  noted 
above  that  the  series  solution  blows  up  at  t/_  =  0,  and  is  not  valid  for  |v_  | 

A  Vtt.  However,  in  the  region  where  the  analytic  solution  is  valid,  that  is  for 
\v_  |  7fc+/fc+  «  0.0035,  the  numerical  and  analytic  curves  are  in  qualitative 
agreement,  in  that  their  central  peaks  have  almost  the  same  width,  and  that 
both  correlation  functions  are  small  when  |v_j  2Av«.  This  confirms  the 
importance  of  the  quantity  Avu  as  the  basic  turbulent  velocity  scale.  We  have 
further  analysed  the  agreement  between  the  two  curves  of  Fig.  (10)  by  plotting 
in  Fig.  (11)  the  width  6vtt  of  the  centred  peak  of  the  correlation  functions,  now 
as  a  function  of  the  parameter  v+,  with  v+  varying  across  almost  the  entire 
range  of  phase  velocities.  The  width  of  the  central  peak  is  defined  as  the 
separation  between  the  two  zeros  at  the  base,  and  is  given  analytically  by  6vtt  = 
4.32 A vti .  The  agreement  remains  rough,  but  reflects  the  correct  dependence 
on  the  resonant  field  amplitudes. 
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While  the  analytic  (7(0,  t/_)  strictly  falls  to  zero  for  |v_|  >  Autt,  the 
numerical  correlation  function  in  Fig.  (10)  retains  sizeable  fluctuations.  We 
suspect  that  this  results  from  the  finite  length  of  the  simulation  (L  =  1024), 
which  limits  the  number  of  uncorrelated  regions  contributing  to  the  spatial 
average  which  defines  the  correlation  function.  For  v+  =  1.745,  the  correlation 
length  is  roughly  lti  ~  90,  and  thus  in  the  neighborhood  of  this  velocity  only 
about  L/ltt  «  10  statistically  independent  regions  exist  in  the  system.  The 
fluctuations  in  Fig.  (10)  are  then  the  noise  which  results  from  doing  statistics 
with  only  about  10  samples.  We  suspect  that  this  noise  could  be  reduced  by 
additional  numerical  processing,  in  which  the  correlation  function,  obtained  as  a 
”  snapshot”  from  Eq.(48),  would  be  further  averaged  in  time,  over  a  period  of  at 
least  a  few  trapping  times  i/^1.  However,  we  have  not  attempted  to  implement 
such  a  scheme. 

While  we  do  not  have  an  analytic  solution  for  (7(0,  v_)  in  any  finite  range 
of  |v_|  A vu  about  v-  =  0,  we  did  find  an  expression  for  the  single  peak 
value  (7(0,0),  Eq.(39).  In  Fig.(12),  we  plot  as  a  function  of  t/+  the  numerical 
and  analytical  values  obtained  for  (7(0,0).  This  graph  indicates  that  Eq.(39) 
is  indeed  qualitatively  correct,  but  that  a  sizeable  quantitative  difference  exists 
between  the  analytic  prediction  of  (7(0,0)  and  the  numerical  result.  We  should 
note  that  the  numerical  value  of  (7(0, 0)  is  relatively  insensitive  to  changes  in 
grid  size  (a  20%  change  results  upon  halving  Ax  or  Av),  a  result  which  gives 
us  confidence  that  the  differences  observed  in  Fig.  (39)  are  not  due  to  spurious 
numerical  effects.  We  suspect  that  the  finite  rate  of  change  of  the  average 
distribution  function,  which  is  not  much  smaller  than  the  growth  rate  (with  the 
estimate  r<f>  «  2000  «4x  7^),  might  be  responsible  for  reducing  the  value 
of  (7(0, 0)  to  under  what  would  be  expected  from  the  uncorrelating  effect  of  7*+ 
alone. 

We  next  consider  the  spatial  structure  of  the  correlation  function,  by  plot¬ 
ting  <7(x_,0),  once  again  for  fixed  t/+.  In  Fig.(13)  this  is  done  for  =  1.745, 
over  the  range  -50  <  x_  <  50.  The  curve  of  <7(x_,0)  displays  the  two  main 
features  predicted  by  the  analytic  theory.  The  first  feature  is  a  relatively  fast 
and  almost  periodic  dependence  on  x_,  with  a  period  we  label  Lpp  (”peak  to 
peak”).  The  second  feature  is  a  broad  envelope,  of  width  we  label  Lm,  which 
modulates  the  fast  dependence  of  the  correlation  function.  In  what  follows, 
Lm  is  defined  graphically,  as  the  length  over  which  the  envelope  decays  to  1/e 
of  its  peak  value,  and  is  estimated  by  linear  extrapolation  from  the  first  three 
central  peaks  of  the  correlation  function.  According  to  the  analytic  results,  we 
should  have  Lpp  =  2tv /fc+  and  «  lu  (eqs.(40)).  The  agreement  is  shown  in 
figs.  (14,15).  It  is  very  good  for  Lpp  and  at  least  qualitatively  correct  for  Lm . 
We  should  stress  however  that  the  result  Lm  «  lit  was  derived  from  very  qual¬ 
itative  considerations;  thus  the  intersection  of  some  parts  of  Fig.  (15)  should  be 
regarded  as  a  coincidence. 

In  Fig. (16) ,  we  plot  the  numerical  C (x_ ,  0)  for  a  smaller  range  of  x_  than  in 
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Fig.(13),  —12  <  x_  <  12,  and  make  a  direct  comparison  with  the  series  solution, 
Eq.(35).  We  note  once  again  that  the  analytic  solution  unfortunately  artificially 
blows  up  at  X-  —  0  (and  at  multiples  of  the  period  2n /k+).  Notwithstanding 
this,  qualitative  agreement  for  values  of  x_  not  too  close  to  these  points  is  good. 


XII.  THE  NUMERICAL  GROWTH  RATES. 

In  this  section  we  consider  the  numerical  results  for  the  growth  rates  which 
were  obtained  in  Simulation  2  (Table  I)  by  the  direct  measurement  of  the  field 
amplitudes.  This  simulation  was  run  with  the  hope  of  observing  a  significant 
enhancement  of  the  growth  rates  over  the  quasilinear  values,  and  in  this  it  differs 
from  Simulation  1  in  two  important  respects:  the  initial  amplitude  spectrum 
\Ek\  was  chosen  uniform  in  fc,  as  opposed  to  \Ek\  ~  1/fc3  previously,  so  as  to 
minimize  nonuniformity  in  trapping  widths,  and  the  system  was  left  to  evolve 
over  a  significantly  longer  time,  that  is  for  about  3.2  e-foldings  of  the  electric 
field  amplitudes.  As  noted  above  (Eq.(45)),  several  e-foldings  should  be  neces¬ 
sary  for  the  equilibration  of  the  nonresonant  harmonics  to  occur  and  hence  for 
the  regime  of  enhanced  growth  rates  to  set  in.  The  initial  linear  growth  rate  was 
7L  =  1.2  10~3  (  obtained  with  a  beam  density  of  Rjb  =  ns /no  =  7.4  10“4)  and 
the  waves  were  initially  excited  in  the  range  of  phase  velocities  0.668  <  v*  <  1.81 
(  77  modes  excited  in  0.552  <  k  <  1.497).  For  the  waves  near  the  central  value 
of  phase  velocity  v*  =  1.0,  we  initially  had  Vtt/~fk  =  5.8,  and  at  the  end  of 
the  run,  i 'tt/lL  =  150,  so  that  as  in  the  previous  simulation,  the  evolution  of 
these  waves  was  in  the  turbulent  trapping  regime  from  the  very  outset  of  the 
instability. 

We  ran  Simulation  2  for  ujpT  =  2875  time  units  (about  450  plasma  peri¬ 
ods),  during  which  the  total  electric  field  energy  increased  about  210-fold.  In 
Fig.(17)  we  show  the  amplitude  spectrum  of  the  forward  waves,  \E^\  at  t  =  0 
and  at  the  final  time  t  =  2875,  that  is  after  about  3.2  e-foldings  of  the  field  am¬ 
plitudes.  The  spectrum  at  t  =  2875  displays  two  distinctive  features.  The  first 
is  that  the  average  spectrum,  which  started  out  as  rigorously  flat,  has  under¬ 
gone  nonuniform  amplification,  with  a  peak  toward  the  lower  values  of  fc,  near 
k  =  0.8.  The  Fourier  amplitudes  for  k  >  0.8  have  not  grown  as  much,  and  there 
has  been  little  growth  at  the  very  edges  of  the  spectrum.  This  nonuniformity 
ran  be  ascribed  to  the  concurrent  flattening  of  the  average  distribution  func¬ 
tion  (Fig.  (18)),  which  is  becoming  severe  towards  the  end  of  the  simulation. 
This  flattening,  when  weighted  with  the  factor  of  v\  in  Eq.(21),  has  a  more 
pronounced  effect  on  the  growth  rate  of  the  waves  with  lower  phase  velocities 
(and  hence  for  the  waves  with  larger  k). 

The  second  distinctive  feature  of  Fig. (17)  is  the  jagged  aspect  of  the  small- 
scale  structure  of  the  spectrum.  There  are  large  and  irregular  amplitude  vari¬ 
ations  from  one  mode  to  the  next,  and  this  results  in  the  spikes  which  can 
be  seen  in  Fig.(17).  We  believe  that  this  irregular  spectrum,  which  was  also 
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observed  by  J.C.Adam  (18),  is  due  to  the  evolution  of  the  system  into  a  set  of 
uncorrelated  subsystems,  each  of  characteristic  length  lu.  Each  Fourier  mode 
can  then  be  written  as  a  sum  of  independent  random  variables ( 24 ),  with,  as  a 
consequence,  the  statistical  independence  of  neighboring  modes. 

To  give  statistical  meaning,  in  terms  of  ensemble  averages,  to  the  mode 
amplitudes  obtained  in  Simulation  2,  we  used  the  convolution  average  defined 
in  Eq.(51)  to  calculate  |£*|2.  We  then  numerically  computed  the  average  growth 
rate  according  to: 


7*  =  “log  |^W|2,  (49) 

which  is  a  definition  of  7 ^  equivalent  to  the  one  assumed  in  Eq.(42).  We  com¬ 
pared  this  to  the  growth-rate  predicted  by  quasilinear  theory,  Eq.(21),  using 
the  numerical,  space  averaged  distribution  function  which  is  obtained  from  the 
simulation: 


7f=  RB\rZ<f'M>„  (50) 

As  a  function  of  velocity,  the  numerical  quantity  <  ff(v)  >t  displays  short-scale 
fluctuations  which  are  probably  the  result  of  the  relatively  short  length  of  the 
simulation,  as  compared  to  the  correlation  length  ( L/lu  «  10).  To  reduce  these 
fluctuations,  we  resorted  to  additional  smoothing  of  <  ff{v)  >,,  by  averaging  its 
value  over  neighboring  points  in  velocity,  through  convolution  with  a  triangular 
window  W(v)  of  the  same  form  as  the  one  used  in  Eq.(51).  While  we  used  the 
same  averaging  function  W  for  <  f(v)  >s  and  |Efc|2,  we  should  note  that  the 
statistical  operations  involved  are  essentially  independent. 

In  Fig.(19),  we  show  the  evolution  of  7jJi(t),  for  k  =  1.0,  vu  =  1.0.  The 
considerable  variation  of  7^  over  the  time  scale  0  <  t  <  2875  is  consistent  with 
the  large  modification  of  <  f(v)  >t  which  can  be  observed  in  Fig. (18),  where  we 
compared  the  initial  and  the  final  average  distribution  functions.  It  is  apparent 
from  this  latter  picture  that  we  are  running  into  a  regime  where  r</>  <  7^_1. 
This  is  a  result  of  the  long  running  time  of  the  simulation,  and  indirectly  of  the 
constraints  placed  on  the  initial  field  amplitudes,  which  cannot  be  infinitesimal, 
but  have  to  be  sufficiently  large  for  the  condition  utt/lfk  >  5  to  be  satisfied  at 
t  =  0. 

Thus,  inspection  of  Fig.  (19)  indicates  that  for  t  >  1000  the  characteristic 
time  r<f>  for  a  change  of  order  unity  to  occur  in  <  /'(t/,£)  >,  is  of  order 
t</>  «  1300,  comparable  to  or  smaller  than  the  e-folding  time  of  the  electric 
field,  which  is  initially  7J“1(0)  =  830,  and  which  decreases  to  7^  «  2000 
towards  the  end  of  the  simulation.  We  may  conclude  that  the  present  simulation 
is  not  an  ideal  verification  of  behavior  in  the  turbulent-trapping  regime,  which 
assumes  . 
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The  evolution  in  time  of  the  ratio  'yk/'1<k  f°r  fc  —  1  is  plotted  in  Fig.(20). 
The  distinctive  feature  of  Fig.  (20)  is  that  7fc/7k'  >  1  throughout  the  evolution 
of  the  system.  A  further  observation  is  that  the  ratio  13  not  constant, 

but  undergoes  a  slow  increase,  reaching  a  value  of  about  1.2  after  2  6  e-foldings 
of  the  electric  field,  at  t  =  2000.  In  the  final  time  interval,  2000  <  t  <  2875, 
Zfk/' y*1  undergoes  a  sharper  rise  reaching  a  maximum  value  of  1.6. 

It  is  tempting  to  explain  the  enhancement  oi^fk/jk  to  values  greater  than  1 
by  the  effects  predicted  by  ALP  and  described  by  the  turbulent-trapping  model. 
However,  we  must  note  that  the  enhancement  of  the  growth  rate  remains  modest 
during  most  of  the  simulation,  with  ^kht  <  l-2  during  the  first  2.6  e-foldings 
of  the  electric  fields.  This  enhancement  is  well  short  of  the  valueof  2.2  predicted 
by  the  turbulent  trapping  model,  Eq.(42).  Furthermore,  when  ~tkhk  increases 
to  a  more  sizeable  value  (from  1.2  to  1.6,  in  the  time  interval  2000  <  t  <  2875), 
it  does  so  precisely  in  a  regime  where  the  trapping  widths  have  become  quite 
large,  that  is  where  neither  quasilinear  theory  nor  the  turbulent  trapping  model 
are  expected  to  be  strictly  valid  in  the  first  place. 

The  correlation  of  large  trapping  widths  with  strong  enhancement  ^  of  the 
growth  rates  is  made  clearer  if  we  look  at  snapshots  of  the  ratio  7^/7*  »  plot¬ 
ted  now  as  a  function  of  fc.  This  is  done  in  figs.(21,22),  for  the  fixed  times 
t  =  2000  and  t  =  2875.  In  the  figures,  we  have  indicated  the  interval  m  k 
corresponding  to  the  turbulent  trapping  width  4Av„,  for  waves  resonant  with 
vk  =  1.0.  Thus,  at  the  end  of  the  simulation  (Fig.(22)),  the  large  enhance¬ 
ment  of  the  growth  rate  occurs  when  the  -rapping  width  is  itself  large,  with 
4Avtt/AvB  ~  0.3.  Furthermore,  this  enhancement  is  accompanied  by  strong 
edge  effects,  so  that  lkhk  is  more  strongly  nonuniform  in  fc  as  well  (similar 
nonuniform  profiles  were  observed  by  J.C.  Adam(18)  ).  This  strong  nonuni¬ 
formity  can  be  ascribed  to  the  fact  that  in  this  regime,  the  waves  derive  their 
energy  from  large  populations  of  trapped  particles,  which  are  trapped  in  only 
a  few  wave  packets  within  the  entire  beam.  The  local  slope  <  f(v)  >,,  which 
defines  the  quasilinear  growth  rate,  is  no  longer  the  salient  physical  parameter 
responsible  for  the  growth  rate  of  the  resonant  waves.  In  this  extreme  regime,  it 
is  not  surprising  that  the  growth  rates  are  uneven,  on  account  of  the  few  wave 
packets  exchanging  energy  with  the  particles,  and  that  they  are  quite  different 
from  the  quasilinear  growth  rates. 

However,  Figs:(20-  22),  have  a  striking  characteristic  which  moves  us  to 
interpret  the  increase  of  growth  rates  as  resulting  from  effects  in  the  spirit  of 
those  predicted  by  ALP.  While  the  nonuniformity  in  the  growth  rate  observed 
in  Fig.  (22)  can  be  ascribed  to  the  large  trapping  widths  of  the  wave  packets, 
this  situation  does  not  by  itself  account  for  the  fact  that  the  growth  rates  are 
greater  than  quasilinear,  everywhere  across  the  spectrum. 

Furthermore,  the  finite  resolution  of  the  grid  severely  limits  the  number 
of  nonresonant  harmonics  which  can  be  generated  by  the  partially  trapped 
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particles.  Thus,  the  maximum  physical  wave  number  k  allowed  on  the  mesh 
is  of  order  tt/2Ax  «  3,  so  that  the  waves  with  k  «  1  may  interact  at  most 
with  3  harmonics.  The  result  of  this  truncation  may  be  a  significantly  smaller 
enhancement  of  the  growth  rate  than  predicted  by  Eq.(42),  which  assumes  that 
an  infinite  cascade  of  harmonics  contibutes  to  the  enhancement.  Thus,  we 
suspect  that  this  may  explain  the  low  enhancements  observed  in  Fig.  (21). 


XIV.  CONCLUSION. 

We  performed  numerical  simulations  to  verify  the  existence  of  a  new  tur¬ 
bulent  regime  in  the  development  of  the  weak  beam-plasma  instability,  and 
more  specifically,  to  validate  the  predictions  of  the  ”  turbulent-trapping”  model 
concerning  correlation  functions  and  enhanced  growth  rates. 

Our  most  definitive  conclusions  regard  the  predictions  of  the  turbulent 
trapping  model  for  structure  in  phase-space.  Both  qualitatively  and  quanti¬ 
tatively,  this  model  provides  good  predictions  for  the  correlation  functions  of 
the  beam  distribution  function.  This  important  result  supports  the  qualitative 
picture  in  which  the  turbulence  is  composed  of  wave  packets  whose  trapping 
domains  barely  overlap  in  phase  space,  with  a  characteristic  width  in  velocity 
A Vft  =  (D/fc)1^3,  and  a  correlation  length  in  space  of  In  =  with 

D^D*1. 

Less  definitive  are  our  conclusions  regarding  the  prediction  of  an  enhance¬ 
ment  of  the  growth  rates,  an  enhancement  which  might  otherwise  be  expected 
on  the  basis  of  the  success  of  the  model  in  predicting  the  correlation  functions. 
The  enhancements  which  are  observed  are  modest  (enhancements  of  1.2-1. 6, 
compared  to  a  prediction  of  2.2),  with  upper  values  occuring  in  a  regime  of 
marginal  validity  for  the  quantitative  predictions  of  the  model.  However,  we 
suspect  that  the  nonideal  aspects  of  the  simulations  are  responsible  for  some  re¬ 
duction  of  the  numerical  growth  rates  which  were  observed.  Thus,  we  conclude 
that  there  is  nonetheless  qualitative  evidence  for  enhanced  growth  rates. 

Future  work  might  bear  on  more  ambitious  simulations,  which  can  over¬ 
come  the  numerical  limitations  of  the  present  ones,  which  suffered  from  a  rela¬ 
tively  coarse  numerical  mesh.  Another  approach  in  the  study  of  the  effects  of 
trapping  on  the  growth  rates  would  be  to  shift  the  emphasis  from  weak  turbu¬ 
lence,  to  the  study  of  what  might  be  called  "moderate”  turbulence  (a  few  wave 
packets  in  the  width  of  the  beam),  a  regime  in  which  the  choice  of  numerical 
parameters  is  less  delicate.  The  insights  gained  from  studying  this  regime  of 
more  vigorous  trubulence  might  then  be  extrapolated  back  to  the  weaker  regime 
assumed  by  the  ALP  model. 


-25- 


Acknowledgements. 

The  authors  wish  to  thank  Dr.  J.C.  Adam  for  suggesting  the  numerical 
approach  taken  in  this  work.  Computation  was  done  at  the  Ecole  Polytechnique, 
Palaiseau,  France,  and  at  the  National  Magnetic  Fusion  Energy  Computing 
Center,  Livermore,  California.  Dr.  K.  Theilhaber  was  funded  on  ONR  contract 
N00014-85-K-0809  and  DOE  contract  FG03-86ER53220. 


-26- 


APPENDIX  A:  ENSEMBLE  AVERAGES  VS.  SPATIAL  AVER¬ 
AGES. 

Statistical  theories  of  turbulence  usually  are  stated  in  terms  of  the  ensem¬ 
ble  averages  of  spectral  quantities.  Thus  one  might  conclude  that  the  proper 
numerical  study  of  turbulence  would  require  performing  a  large  number  of  nu¬ 
merical  simulations,  each  with  different  initial  conditions,  and  averaging  the 
final  results.  In  fact,  if  the  turbulence  studied  is  homogeneous  in  space,  in 
principle  all  the  relevant  information  about  ensemble  averages  can  be  extracted 
from  any  one  realization  in  the  ensemble,  that  is  from  a  single  numerical  sim¬ 
ulation,  provided  that  the  system  is  long  enough.  The  criterion  for  sufficient 
length  is  that  the  length  L  of  the  system  be  much  longer  than  the  correlation 
length  lc  of  the  quantity  investigated  (see  (25)  for  equivalent  discussions  of  data 
processing,  for  power  spectra  in  the  time  and  frequency  domain). 

For  instance  consider  the  squared  Fourier  amplitude  of  the  electric  field, 
which  will  be  the  main  spectral  quantity  of  interest.  The  problem  is  to  approx¬ 
imate  its  ensemble  average,  which  we  denote  by  <  \Ek\ 2  >,  from  the  results  of 
a  single  simulation.  The  prescription  is  as  follows:  we  first  obtain  the  ”raw” 
squared  Fourier  spectrum  |Efc|2,  directly  from  the  numerical  E(x)  resulting 
from  a  single  computer  run.  In  general,  this  quantity  will  be  a  rather  chaotic 
function  of  k.  We  then  perform  a  convolution  to  define  the  smoothed  average: 


\Ek\2  =  Y,W(K~k)\EK\2,  (51) 

K 

where  W{K)  is  a  window  function  such  that  W(K)  =  1.  In  what  follows, 
we  choose  W(K)  to  be  triangular  in  shape,  with  W(K)  =  (Afc/fcv,)(l  —  |If|  /fc^) 
and  W(K)  —  0  for  \K\  >  k^.  Provided  that  the  width  of  W(K )  is  small 
compared  to  the  total  spectral  width  (that  is,  the  inverse  correlation  length), 
but  large  compared  to  the  mode  spacing  (&**  A k  =  2i r/L),  Eq.(51)  will 

provide  an  approximation  to  the  ensemble  average,  \Ek\2  ~<  \Ek\2  >. 

In  dealing  with  spatial  quantities,  such  as  the  spatial  correlation  functions 
or  the  average  distribution  function,  it  is  legitimate  to  replace  the  ensemble 
average  by  a  sliding  spatial  average  over  the  entire  length  of  the  system,  pro¬ 
vided  once  again  that  the  system  is  long  compared  to  the  correlation  length.  In 
what  follows  we  adopt  the  following  conventions:  the  unadorned  angle  brackets 
<  . . .  >  refer  to  ensemble  averages,  the  addition  of  a  subscript  ”s”  denotes  a 
spatial  average,  <  ...  >,,  and  the  overbar  notation  will  be  reserved  for  the 
spectral  quantities,  to  denote  the  filtering  operation  of  Eq.(51). 
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APPENDIX  B:  CONDITIONS  OF  VALIDITY c 

In  this  Appendix  we  summarize  the  several  conditions  which  define  the 
turbulent  trapping  regime.  In  addition  to  those  stated  in  the  main  text,  a 
condition  on  the  validity  of  the  derivation  of  (42,43)  is  that  dispersion  within  a 
turbulent  wave  packet  is  small  over  a  time  scale  of  1/7*.  This  is  expressed  by 
the  condition  (4): 


(Puk 

dk 2 


( 


Vtt 

\Vk-Vgk  I 


) 


2 


<7fc, 


(52) 


where  v*.  is  the  phase  velocity  and  vgk  the  group  velocity  of  the  waves.  The 
cold  bulk  plasma  is  theoretically  dispersionless,  but  in  fact  has  the  numerical 
dispersion  relation  u>k  =  1  —  fc2 Ax2/ 8  valid  for  kAx  -C  1.  With  v *  «  1/fc, 
=  \duk/dk\  Vk ,  Eq.(52)  reduces  to: 


_  1  fc2AxV< 

=  7 - 1, 

4  7^ 


(53) 


This  constraint  is  necessary,  in  addition  to  those  already  outlined,  for  a  valid 
test  of  the  turbulent  trapping  theory. 

To  summarize  all  constraints,  we  have  (eqs.(24),(41),(44)  and  (53),  and  the 
paragraph  following  Eq.(44)): 


Tc  <  Vu1  <7fc1  £T</>» 

(54) 

4Avtt  -C  A vb, 

(55) 

vt  t  >  5  7k, 

(56) 

_  1  k2Ax2tjt  „ 

*disp  —  .  ^  1) 

4  7k 

(57) 

TV  «  3-5, 

(58) 

where  iV  is  the  number  of  e-foldings  during  the  evolution  of  the  instability. 
Eqs. (54,55)  are  general  conditions  for  weak  turbulence,  and  are  required  for 
quasilinear  theory  as  well.  Only  eqs.  (56,57,58)  specifically  define  the  turbulent 
trapping  regime. 

Eqs. (55),  (56)  and  (58)  are  the  most  stringent  conditions  imposed  on  the 
simulations.  For  instance,  suppose  that  we  want  to  follow  the  evolution  of  the 

system  through  N  e-foldings  of  the  field  amplitudes,  while  remaining  in  the 
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turbulent  trapping  regime  throughout  the  computer  run.  For  a  fixed  beam 
width,  Eq.(55)  limits  the  amplitude  of  the  final  electric  field.  This  implies  that 
the  initial  field  amplitudes  have  to  be  scaled-down  by  a  factor  of  order  e~N  from 
a  rigid,  fixed  upper  limit.  Thus  the  initial  value  of  uit  «  \E\2^3  «  e-2JV/3  is 
small.  In  turn,  Eq.(56)  implies  that  the  initial  growth  rate  must  be  sufficiently 
small  as  well,  with  the  same  scaling  7  ~  vuf  5  «  e~2N^3.  Thus  the  total  run 
time  T  of  the  simulation,  which, according  to  Eq.(58),  has  to  undergo  at  least  N 
e-foldings,  will  scale  a sTk  JVe2^3.  We  conclude  that  the  computation  time 
increases  very  rapidly  with  the  number  of  e-foldings  we  wish  to  observe.  The 
value  N  =  4  represented  a  practical  limit  to  our  computing  resources. 

APPENDIX  C:  INITIALIZATIONS. 

The  initial  beam  distribution  function  was  designed  for  a  constant  linear 
growth  rate  over  almost  all  of  the  unstable,  positive-slope  velocity  range,  by 
piecing  together  the  simple  functions  ( 19 ): 


A(v  - 1/1)2  +  B( v  —  vx)\  vi  <  v  <  va, 
fo(v)  =  Ci  -I-  C2(l  -  va/v ),  va  <  v  <  vb ,  (59) 

D(v  —  V2)2  -f  E(v  —  V2)3,  vb  <  v  <  V2 

In  Eq.(59),  the  expressions  in  the  velocity  ranges  (vi,  va)  and  (uj,,  V2)  are  transi¬ 
tion  functions,  with  the  coefficients  chosen  so  that  /o(v)  and  fo(v)  are  continu¬ 
ous  everywhere.  The  coefficients  are  further  chosen  to  satisfy  the  normalization 
/  fo{v)dv  —  1.  The  form  (59)  insures  that  in  the  interval  va  <  v  <  vby  the 
growth  rate,  given  by  Eq.(21),  is  constant. 

A  concern  in  choosing  the  coefficients  e*  in  Eq.(46)  is  to  insure  that  the 
resonance  overlap  between  adjacent  modes  is  satisfied  within  the  entire  excited 
spectrum,  and  this  from  the  very  beginning  of  the  simulation.  Considerations  of 
trapping  widths  axe  important  because  in  the  absence  of  resonance  overlap  there 
is  no  diffusion.  In  other  words,  the  effect  of  the  discreteness  of  the  spectrum 
is  attenuated  only  provided  the  modes  overlap.  In  fact  resonance  overlap  was 
satisfied  in  all  of  our  simulations. 
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TABLE  1 :  Simulation  Parameters . 
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FIG.  1.  Sketch  of  the  bulk  and  beam  distribution  functions. 
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FIG.  2.  (a)  Plot  of  Zf(0,u)  showing  the  velocity  dependence  of  the  analytic  correlation 

function;  u  =  v_  /2l^2Avtt.  (b)  Plot  of  H(£, 0)  showing  the  spatial  dependence  of  the 
analytic  correlation  function;  £  =  k+x^. 
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FIG.  2.(c)  Contour  plots  of  the  analytic  correlation  function  H{^u). 
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FIG.  3.  (a)  Initial  beam  distribution  function  /o(t>)  for  Simulation  1  and  phase  velocities 

of  the  Fourier  components  excited  at  t  =  0.  (b)  Initial  linear  growth  rate  for  Simulation  1. 
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FIG.  5.  Final  space-averaged  distribution  function  for  Simulation  1;  the  dashed  line  is 
the  initial  space-averaged  distribution  function,  fo{v). 

FIG.  6.  A  cross-section  of  the  final  distribution  function  in  simulation  1,  taken  at  x  =  512, 
showing  fluctuations  in  /;  the  dashed  line  is  the  final  space- averaged  distribution  function. 


FIG.  7.  Contour  plot  of  the  final  beam  distribution  function,  at  t  =  500,  for  Simulation 
1.  The  total  system  length  is  1024. 
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FIG.  8.  Perspective  view  of  the  beam  distribution  function  in  Area  1  of  Fig. (7).  Note 
the  plateaus  which  have  formed  in  the  distribution  function. 
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FIG.  9  Perspective  view  of  the  beam  distribution  function  in  Area  2  of  Fig. (7), 
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FIG.  10.  Comparison  of  the  numerical  and  analytic  C(0,t>_)  for  Simulation  1,  for  the 
mean  velocity  r+  =  1.745,  at  t  —  500. 

FIG.  11.  Comparison  of  the  width  in  t>_  of  the  function  C(0,  v_)  (the  separation  of  zeros 
at  the  base),  plotted  as  a  function  of  the  mean  velocity,  for  Simulation  1,  t  =  500.  The 
analytic  results  are  obtained  with  Svu  =  4.32Autf. 
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FIG.  12.  Peak  values  C(0,0)  obtained  in  Simulation  1,  t  =  500. 
FIG.  13.  Plot  of  C(x_,0)  for  =  1.745,  Simulation  1,  t  =  500. 


FIG.  14  Comparison  of  the  periodicity  length  Lpp  in  C(x_,0)  from  Fig. (13).  The  ana¬ 
lytic  results  are  obtained  with  Lpp  =  2n /fc+, 

FIG.  15.  Comparison  of  numerical  and  theoretical  spatial  decay  lengths  for  C(x_,0)  from 
Fig. (13),  The  analytic  result  is  obtained  with  ltt  =  l/k+ZD1/2. 
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FIG.  16.  Plot  of  C(x_,0)  for  t>+  =  1.745,  Simulation  1.  This  is  a  blow-up  of  Fig. (13). 
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FIG.  17.  The  amplitude  spectrum  \E£\  in  Simulation  2,  a)  at  t  =  0  and  b)  at 
t  =  2875.  The  dashed  line  denotes  the  Fourier  amplitudes  at  t  =  0. 

FIG.  18.  The  initial  and  final  average  distribution  function,  as  obtained  m  Si nmlfttioD  2. 


FIG.  19.  Evolution  of  the  quasilinear  growth  rate  of  the  Fourier  component 
fc  =  1.0,  Vk  —  1.0,  in  Simulation  2,  as  calculated  from  eq.(50). 


FIG  20.  Evolution  of  the  ratio  / 7^  in  Simulation  2,  for  A;  =  1.0.  The  two 
vertical  numbers  indicate  the  approximate  number  of  e-foldings  undergone  by 
the  field  at  a  given  time. 


FIG.  21.  The  ratio  of  the  growth  rates  in  Simulation  2,  Ik/lk!’  plotted  as  a 
function  of  k  for  t  =  2000.  The  interval  indicates  the  range  in  k  covered  by 
the  trapping  width  of  waves  near  k  =  1.0(the  edges  of  the  spectrum  have  been 
excluded  by  the  averaging  process  which  defines  7fc). 

FIG.  22  Same  as  Fig. (5),  but  for  t  =  2875. 
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